Variable types

Before diving into regression types, let’s briefly review variable types:

• Quantitative variables (numerical values) can be:
  • Discrete: countable values (e.g., number of wind turbines in district)
  • Continuous: any value within a range (e.g., income)
• Qualitative (categorical) variables (non-numerical) can be:
  • Nominal: categories without order (e.g., env_concern)
  • Ordinal: categories with a meaningful order (e.g., educ)

A binary variable is a special case of nominal with only two outcomes (e.g., yes/no).

Working with categorical data with {forcats}

In R, factors are used to represent categorical variables. A factor is a compound object made up of two components:

• A vector storing the actual category values
• A set of levels defining the possible categories (they can also take ordering)

# creating a factor
countries_of_nationality <- factor(
  c("Italy", "France", "Germany", "Germany", "Austria")) 

countries_of_nationality

## [1] Italy   France  Germany Germany Austria
## Levels: Austria France Germany Italy

# creating a factor with ordered levels (see what happens when you remove ordered = TRUE)
sizes <- factor(
  c('sm', 'md', 'lg', 'sm', 'md'),
  levels = c('sm', 'md', 'lg'),
  ordered = TRUE)

sizes

## [1] sm md lg sm md
## Levels: sm < md < lg

While R has a native way to work with factors, the {forcats} package from the tidyverse provides a useful and flexible alternative. Write more pitch for it. {forcats} simplifies working with factors, especially when dealing with:

• Reordering levels: Functions like fct_relevel() allow you to reorder factor levels easily, useful for controlling the order of categorical variables in analysis or plotting.
• Recoding levels: fct_recode() allows you to recode levels, merging or renaming categories without needing to manually manipulate the data.
• Handling missing values: {forcats} provides functionality to explicitly handle missing values in factors using fct_explicit_na().
• Level manipulation: With functions like fct_drop(), it is simple to drop unused levels, which can be crucial when cleaning data.
• Summary and visualization: {forcats} integrates seamlessly with ggplot2 and dplyr to summarize and visualize categorical data efficiently.

Some of the most useful functions are:

• forcats::fct_reorder(): Reordering a factor by another variable.
• forcats::fct_infreq(): Reordering a factor by the frequency of values.
• forcats::fct_relevel(): Changing the order of a factor by hand.
• forcats::fct_lump(): Collapsing the least/most frequent values of a factor into “other”.

We will use some of them later today.

You can take a look at the chapter on Factors in the R for Data Science book and the {forcats} vignette here.

categorical_variables <- c("support_wind_turbines","proximity_to_proposed",
                           "env_concern","trust_in_government",
                           "left_right_ideology", "educ")

wind_factors_df <- wind_df %>% 
  dplyr::mutate_at(dplyr::vars(categorical_variables), as.factor) %>% # turn into factors
  dplyr::mutate(
    educ = forcats::fct_recode(
      educ,
      "High school or less" = "1",
      "University degree" = "2",
      "Postgraduate degree" = "3"
    ),
    env_concern = forcats::fct_recode(
      env_concern,
      "Low" = "1",
      "Medium" = "2",
      "High" = "3"
    ),
    income_1000 = income/1000) #income in 1000s

#wind_factors_df$env_concern

Linear regression

• Simple linear regression: one independent variable, continuous dependent variable.
• Multiple linear regression: two or more independent variables, continuous dependent variable.

Logistic regression

• Binary logistic regression: binary dependent variable (two outcomes).
• Multinomial logistic regression: nominal dependent variable with 3+ unordered categories.
• Ordinal logistic regression: ordinal dependent variable with 3+ ordered categories.

Logistic and probit regressions, which you saw during this week’s lecture, are both part of Generalized Linear Models (GLM), which extend linear models to handle non-continuous outcomes and non-normal residuals. In the labs we will focus on logistic regression, but the general syntax to set up the models and logic of extracting predictive margins with {marginaleffects} can be applied to probit models.

Logistic regression with one independent variable

Say we want to examine the relationship between income and the likelihood of supporting wind turbine projects.

A binary logistic regression is appropriate when the dependent variable is categorical with two levels (e.g., support vs. no support), and the independent variable(s) can be either categorical or continuous.

In this case, the dependent variable is whether a resident supports the wind turbines (coded as 1 for support and 0 for no support), and the independent variable is income (a continuous variable).

In R, logistic regression can be done using the glm() function, specifying with the argument family = "binomial". The syntax is similar to other formulas we have used in the past:

glm(outcome_of_interest ~ predictor,
family = binomial,
data = data_frame)

income_model <- 
  glm(support_wind_turbines ~ income_1000,
      family = binomial,
      data = wind_factors_df
  )

modelsummary::modelsummary(income_model,
                           stars = c('*' = .1, '**' = .05, '***' = .01),
                           statistic = 'conf.int')

exp(0.039)

## [1] 1.03977

predict(income_model,
        newdata = data.frame(income_1000 = 70),
        type = "response"
)

##         1 
## 0.8805791

pred_income <- predict(income_model,
                       newdata = data.frame(income_1000 = 70),
                       type = "response",
                       se = T
)

pred_income$fit

##         1 
## 0.8805791

## creating CIs
lower <- pred_income$fit - (qnorm(0.975) * pred_income$se.fit)
upper <- pred_income$fit + (qnorm(0.975) * pred_income$se.fit)
c(lower, upper)

##         1         1 
## 0.8747801 0.8863782

marginaleffects::plot_predictions(income_model, condition = "income_1000") 

marginaleffects::plot_predictions(income_model, condition = "income_1000") +
  geom_point(x = 70, y = 0.8805791, color = "#cc0065", size = 2) + 
  geom_vline(xintercept = 70, alpha = 0.5, color = "#cc0065", linetype = "dashed") + # this is our prediction for 70
  theme_minimal() +
  labs(x = "Income (thousands)",
       y = "Predicted probability of support\nfor wind turbines")

Logistic regression with multiple independent variables

The interpretation of the coefficients in multivariable logistic regression is similar to the interpretation in univariable regression, except that this time it estimates the multiplicative change in the odds in favor of \(Y = 1\) when \(X\) increases by 1 unit, while the other independent variables remain unchanged.

For this illustration, suppose we would like to estimate the relationship between supporting wind turbine projects and all variables present in the data frame:

multivariate_model <- glm(support_wind_turbines ~ proximity_to_proposed + env_concern +
                            trust_in_government + left_right_ideology + age + income_1000 + educ,
                          family = binomial,
                          data = wind_factors_df)

modelsummary::modelsummary(multivariate_model,
                           stars = c('*' = .1, '**' = .05, '***' = .01),
                           statistic = 'conf.int')

Like previous example with one predictor, it is easier to interpret these relationships through OR. But this time, we also print the 95% CI of the OR in addition to the OR (rounded to 3 decimals) so that we can easily see which ones are significantly different from 1:

round(exp(cbind(OR = coef(multivariate_model), confint(multivariate_model))), 3) %>% 
  knitr::kable() %>% 
  kableExtra::kable_styling()

From the ORs and their 95% CIs computed above, we conclude that:

• Proximity to the proposed site: individuals who live close to the proposed site have odds of supporting the project that are 0.46 times those of individuals who do not live nearby, all else being equal. This indicates a negative association between proximity and support.

• Environmental concern: individuals with medium environmental concern have odds of support that are 3.31 times those with low concern, while individuals with high environmental concern have odds 10.24 times higher than those with low concern, all else being equal. These strong and significant associations suggest that higher environmental concern is positively associated with support for the project.

• Trust in government: individuals who trust the government have odds of supporting the project that are 1.69 times higher than those who do not, all else being equal.

• Left-right ideology: individuals who lean right on the ideological spectrum have odds of support that are 0.52 times those who lean left, indicating a negative relationship between right-leaning ideology and support for the project.

• Age: for each additional year of age, the odds of supporting the project are multiplied by 0.98, suggesting that older individuals are slightly less likely to support the project, all else being equal.

• Income: for each additional 1,000 Burovian dollars of income, the odds of supporting the project increase by a factor of 1.05. This positive and significant association indicates that wealthier individuals are more likely to support the project.

• Education: individuals with a university degree have odds of support that are 1.34 times those of individuals with a lower level of education, while those with a postgraduate degree have odds that are 1.89 times higher, all else being equal.

• Intercept: we refrain from interpreting the intercept, as it represents the odds of support for a hypothetical individual with all predictors set to zero or to the reference category, which is not substantively meaningful in this context.

Interactions

In the previous sections, we omitted potential interaction effects.

As we learned last week, an interaction occurs when the relationship between an independent variable and the dependent variable depends on the value taken by another independent variable. In other words, if the effect of one predictor on the outcome changes depending on the level of another predictor, then there is evidence of an interaction.

Based on what we know from Burovia, we hypothesize that our expectations of proximity to a proposed development on support for the development may vary depending on a respondent’s level of environmental concern.

interaction_model <- glm(support_wind_turbines ~ proximity_to_proposed * env_concern +
                            trust_in_government + left_right_ideology + age + income_1000 + educ,
                          family = binomial,
                          data = wind_factors_df)

modelsummary::modelsummary(interaction_model,
                           stars = c('*' = .1, '**' = .05, '***' = .01),
                           statistic = 'conf.int')

marginaleffects::avg_predictions(interaction_model, by = c("env_concern", "proximity_to_proposed"))

## 
##  env_concern proximity_to_proposed Estimate Std. Error     z Pr(>|z|)   S 2.5 %
##       Low                        0    0.667    0.00535 124.6   <0.001 Inf 0.656
##       Low                        1    0.529    0.00860  61.6   <0.001 Inf 0.513
##       Medium                     0    0.823    0.00525 156.8   <0.001 Inf 0.812
##       Medium                     1    0.734    0.00892  82.2   <0.001 Inf 0.716
##       High                       0    0.926    0.00381 243.0   <0.001 Inf 0.918
##       High                       1    0.866    0.00735 117.8   <0.001 Inf 0.852
##  97.5 %
##   0.677
##   0.546
##   0.833
##   0.751
##   0.933
##   0.880
## 
## Type:  response 
## Columns: env_concern, proximity_to_proposed, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

marginaleffects::plot_predictions(interaction_model, by = c("env_concern", "proximity_to_proposed")) +
  labs(x = "Environmental concern level",
       y = "Predicted probability of support\nfor wind turbines",
       color = "Proximity to proposed site") +
  scale_color_manual(labels = c("0" = "Does not live nearby", 
                                  "1" = "Lives nearby"),
                       values = c("#5C9184","#550055")) +
  theme_minimal() +
  theme(legend.position = "bottom")